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(73 Abstract 

A first-principle multiscale modeling approach is presented, which is derived from 

^ the solution of the Ornstein-Zernikc equation for the coarse-grained representation of 

1^ polymer liquids. The approach is analytical, and for this reason is transferable. It is 

P3 here applied to determine the structure of several polymeric systems, which have differ- 

2 ent parameter values, such as molecular length, monomeric structure, local flexibility, 

I— —I and thermodynamic conditions. When the pair distribution function obtained from 

^ this procedure is compared with the results from a full atomistic simulation, it shows 

J> quantitative agreement. Moreover, the multiscale procedure accurately captures both 

O large and local scale properties while remaining computationally advantageous. 

in 
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o 

g 1 Introduction 

Understanding the structure of complex fluids hinges on our ability of achieving a detailed 
^ picture of the system on an extended range of lengthscales. Many important phenomena, 

^ which are tightly dependent on the local chemical structure, occur on a large scale. Ther- 

modynamic properties, which can be considered as characterized by a macroscopic (infinite) 
lengthscale change as the local, monomeric structure is modified. For example compressibil- 
ity or phase diagrams of a macromolecular liquid depend on the type of polymer involved. 
As the physics appearing on different lengthscales is related, an effort has to be made to 
formally correlate information on different lengthscales.^'^ 

Macromolecular liquids are complex, given the plethora of characteristic lengthscales 
that emerges from their investigation. Two fundamental, molecular lengthscales are relevant 
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in the study of any macromolecular liquid, namely the monomer lengthscale, I, which is 

local, and the molecular radius-of-gyration, Rg = l^J N/Q, where is the molecular degree 
of polymerization, which roughly represents the overall size of the polymer. Because Rg 
increases linearly with the number of monomers in the chain, those two lengthscales can 
be separated by many orders of magnitude in real macromolecular systems of practical 
application. 

The local scale is an effective length in which complicated local information is embed- 
ded. For example, in the simplest and most conventional description of the polymer local 
structure, inter- and intra-molecular forces balance each other leading to an unperturbed 
polymer chain statistics at any lengthscale.^ However, simulations show that on the local 
scale entropic effects and chain connectivity produce a segmental correlation hole that man- 
ifests itself as local, intramolecular clustered, structures of monomers.^ Because the liquid 
is almost incompressible and monomers from the same chain tend to cluster together due to 
connectivity, the ideality hypothesis does not hold on the local scale. 

A second effect important on the local scale is the presence of rotational energy bar- 
riers.^' ^ To describe liquids of semiflexible polymers, where the local dynamics is affected 
by conformational transitions, the bond length is modified into an effective segment length 
that includes the effect of hindered rotations. This lengthscale is commonly referred to as 
the polymer persistence length.^ 

Computer simulations are useful in providing detailed information about the macro- 
molecular systems under study. Correlation functions can be readily calculated from simu- 
lation trajectories and compared with experiments. However, simulations are limited in the 
range of lengthscales they can cover, as the statistics is good up to the lengthscale defined 
by half the size of the simulation box. Increasing the box size, increases the number of 
interacting particles and the number of pair potentials that need to be calculated at any 
timestep, dramatically slowing down the simulation. 

For a system that is ergodic, where statistical averages over an ensemble of particles is 
equivalent to the average of one system over many computational time steps, one could think 
of limiting the number of particles simulated, while still being able to recover enough statistics 
by extending the length of the simulation run. Unfortunately, the problem of error build-up 
at each step during the simulation run makes calculated trajectories rapidly unreliable. This 
is due to the fact that trajectories are derived by numerical solutions of differential equations. 
These solutions are chaotic in nature, and small differences due to numerical uncertainty in 
the initial conditions produce large fluctuations in the output trajectory, which increase with 
time. The extent to which simulation trajectories are sensitive to small initial errors depends 
on the value of the Lyapunov exponent. Shadowing trajectories, which well represent the 
real ones, only last for a limited number of computational steps. ^'^ 

Other methods have been studied to improve the maximum number of computational 
steps that can be reliably simulated, including hybrid methods, e.g. implicit solvent methods, 
that conserve the atomistic description locally while assuming a meanfield description outside 
a local volume. However, these methods remain under discussion because of the errors that 
they can include in the procedure.^ 

In a nutshell, if the system is described at high resolution, a large number of instan- 
taneous forces have to be calculated at each computational time step, which slows down 
the computational time. Computational time steps have to be short to avoid large initial 
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errors, but because trajectories become increasingly unreliable with the number of steps, this 
limits the longest timescale that can be reached in the simulation, and so the quality of the 
simulation output/'^ 

A general strategy useful for increasing the timescale and the number of particles that 
can be simulated is to adopt a coarse-grained description of the system. The simulation 
of a coarse-grained system can include a significantly larger number of particles and can 
extend to a much longer time scale than simulations of the same system described with 
atomistic resolution. This is because internal degrees of freedom and local energy barriers 
are averaged out, allowing a larger elementary time step than for atomistic simulations. 
Moreover, in the coarse-grained description each molecule contains a number of atoms that 
can be considerably reduced with respect to the full atomistic descriptions, so that the 
number of pair interactions that have to be calculated, at each computational step between 
each pair of molecules, is reduced, allowing for the simulation of larger systems. 

The tradeoff of coarse-grained methods is that the local information is completely lost, 
and also that because the local energy barriers are averaged out and internal degrees of 
freedom are eliminated, the simulated dynamics is orders of magnitude faster than the real 
one, and there are not easy methods to recover the correct dynamics.^'' 

Modeling of coarse-grained systems have recently received considerable attention^^ as 
simulations aspire to describe more and more complicated structures relevant for engineering 
applications as well as biophysical systems. ^^"^^ The key question in developing coarse- 
grained molecular models is the definition of the effective potential acting between units. The 
general strategy pursued to solve this problem is a bottom- up approach. First, atomistic 
simulations of the molecular liquid are performed. Once the coarse-grained structure is 
defined by lumping together adjacent atoms in the molecule, an effective numerical coarse- 
grained potential is obtained, which is parametrized by optimizing the agreement against 
properties measured in the atomistic simulation. For example, the parameters in the effective 
potential can be defined by the optimization of the mean-force potential, matching the forces, 
or by optimizing the free parameters against known physical properties. The potential, so 
derived, is used as an input to a new simulation of the coarse-grained system. The properties 
derived from this simulation are checked for consistency against physical properties, e.g. pair 
distribution functions, obtained from the atomistic simulation, and the process is repeated 
until a fully self-consistent description is achieved. The test quantity that is optimized is in 
most cases the pair distribution function, or analogously the total correlation function, as 
from this quantity all the properties of a liquid can be directly calculated. 

This numerical procedure has the limitation that the resulting optimized effective po- 
tentials are specific to the chemical structure considered and of the thermodynamics condi- 
tions in which the initial atomistic simulations were performed. In principle those potentials 
cannot be applied to describe the system in a different thermodynamic condition, which 
means that the potential is not transferable. 

Several coarse-graining methods have been presented in the literature. An extremely 
successful approach has been the United Atom (UA) description where each moiety of type 
C, CH, CH2, or is represented as an effective unit.^^ Effective potentials between 

united atoms have been derived and parametrized to reproduce the properties of the system 
described at the atomic level. The UA description is useful because polyolcfins, which 
constitute a great portion of the macromolecular systems of interest, can be fully represented 
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as a collection of these sites. The UA description has been extensively tested and has been 
proven to be extremely reliable. ^^"^^ 

Combined with new techniques that speed up the equilibration step, the UA description 
allows for MD simulations of polymer melts with long entangled chains. 

While UA simulations have proved successful in reproducing the structural and dynam- 
ical properties of polymers, there is an increasing need to extend simulations to even larger 
scales to capture other phenomena, such as cooperative self-assembly of supermolecular struc- 
tures, dynamics of systems approaching their second order phase transitions, simulations of 
bulk properties, to name a few. Recently, we have developed an original coarse-grained 
method, which is analytical and represents polymers on the lengthscale of Rg as interpene- 
trating soft colloidal particles. ^^"^^ A more refined coarse-grained description has also been 
derived, where each macromolecule is represented as an elastic dumb-bell of interpenerating 
spheres. The effective potential in these descriptions is derived by solving the Ornstein- 
Zernike equation in the monomer and effective unit sites, where the center-of-mass (com) 
of the effective unit is assumed to be an auxiliary site, i.e. there is no direct correlation 
between monomer and com. Relaxing this approximation does not change the outcome of 
the coarse-graining procedure. Because the monomer intramolecular distribution follows 
a Gaussian statistics, the Ornstein-Zernike equation can be solved analytically to produce 
the pair distribution function for the effective units. From the latter, the effective coarse- 
grained potential is calculated by applying the hypernetted-chain closure approximation.^^ 
The advantage of our coarse-grained description is that the obtained potential is directly 
transferable to systems in different thermodynamic conditions and to polymers with differ- 
ent degree of polymerization or different chemical structure. ^^"^^ Comparison against UA 
simulations of several polymeric liquids and mixtures showed quantitative agreement between 
the total correlation functions in the two representations. 

Because the trade-off of coarse-graining procedures is that the information on the local 
length scale is lost, techniques have been developed to merge the coarse-grained descrip- 
tion with local information. In the so-called "multiscale modeling" procedures a hierarchy 
of simulations are performed on the system coarse-grained at different characteristic length- 
scales, and then the information obtained from those simulations is combined into one overall 
description, which includes all the lengthscales of interest. 

The method that we are presenting here is a multiscale modeling procedure for the 
simplest macromolecular liquids, i.e. of homopolymers. Here, the lengthscales of interest 
are the monomer and the radius-of-gyration, so that two simulations have to be combined, 
one at the monomer level of description, and one where polymers are represented at the 
Rg level, or mesoscopic simulation (MS-MD). We perform molecular dynamic simulations 
of the liquid of polymers described as soft colloidal particles in a MS-MD simulation and 
we combined those results with the outcome of a united atom molecular dynamic (UA-MD) 
simulation for the local description, through our multiscale modeling procedure. We test 
our procedure for a number of systems, which include polyethylene melts with different 
degrees of polymerization, as well as liquids of macromolecules with different monomeric 
structures. The procedure is optimized and shows quantitative agreement between the total 
correlation functions calculated by us, and the ones obtained in a test UA-MD simulation. 
The advantage of our procedure, with respect to running a full UA-MD simulation of the 
same system, is due to the gain in computational time. 
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Another advantage of our multiscale modeling procedure, with respect to other mul- 
tiscale approaches, is that it is formally compatible with the first principle formalism used 
to coarse-grain the polymeric structure as both rely on the Ornstein-Zernike equation which 
includes auxiliary sites. Because it is analytical, the coarse-graining procedure can be used 
not only to remove internal degrees of freedom, but also in the opposite way to reintroduce 
"a posteriori" those internal degrees of freedom after the MS-MD is completed. It is shown 
that this approach reproduces pair correlation functions at a high computational efficiency, 
providing a method of extending simulations to large length scales of interest. 

2 Coarse-graining of a polymer at the radius-of-gyration 
lengthscale 

Our coarse-grain model, briefly revisited here, maps each macromolecule into an interacting 
soft-colloidal particle. We developed this model in a series of papers, where we reported 
coarse graining of polymer melts with different architectures as well as coarse graining of 
mixtures of those polymers.^' ^^"^^'^^ 

The first modeling of polymers as soft colloidal particles dates back to Flory and Krig- 
baum,^^ where repulsive monomer-monomer interactions are integrated over the Gaussian 
chain distribution. The resulting effective potential at contact is repulsive and becomes 
stronger with increasing chain length or polymer density, in disagreement with simulations, 
scaling theories and renormalization group calculations. Correct scaling behavior was ob- 
tained by Hansen and coworkers who derived pair distribution functions for polymers in 
dilute or semi-dilute solution, starting from first principles liquid state theory.^^ A more 
phenomenological form of this approach was implemented early on to describe polymer melts 
and mixtures by Dautenhahn and Hall,^® and by Murat and Kremer.^^ 

The choice of the monomer and radius-of-gjnration as the two length scales at which 
the polymer liquid is coarse-grained in our multiscale procedure is motivated by the fact that 
the monomer pair distribution function, and the structure of the polymer liquid, are fully 
determined when these two length scales are known, as shown below. 

Our coarse-graining procedure starts from Curro and Schweizer's polymer reference 
interaction site model (PRISM) Ornstein-Zernike equation, which for a homopolymer melt 
relates the monomer total pair correlation function, h{r), to the direct correlation function, 
c(r), in reciprocal space^*^''^^ 

h{k) = u{k)c{k)[u{k) + ph{k)] , (1) 

where i2}{k) is the intramolecular structure factor, and p is the monomer site number density. 
The procedure is then extended to include the center-of-mass of the coarse-grained structures 
as auxiliary sites, using the formalism of Krakoviack et al.,^"^ which starts from a matrix form 
of Eq.Q. The density of auxiliary sites is given by the molecular number density, pch = p/N . 
The solution of Eq.([T]) is obtained under the assumptions that the direct correlation functions, 
c(/c), only depend on interactions between real sites, while no interaction is assumed between 
auxiliary sites. Relaxing this condition does not modify the zeroth order result, ^'^ which yields 
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the equation 



UJ'' 



(2) 



where h'^"^{k) is the total pair monomer-monomer correlation function, uj'^'^{k) is the in- 
trachain correlation function of monomers about the com, a;'"™(/c) is the intramolecular 
monomer-monomer correlation function, and h'^^{k) is the mesoscale com-com total correla- 
tion function. Eq.(|2| is the basis for our coarse graining approach since it links the com and 
the monomer site correlation functions. The com representation corresponds to a mesoscale, 
coarse-grained description of the polymeric liquid, where each chain is a soft colloidal par- 
ticle centered on the position of the polymer com. On the other hand, the monomer site 
description is related to the united atom representation of the polymer chain, where each 
site represents an united atom. 

In order to obtain an expression for the coarse-grained description, we adopt for the 



monomer correlation function, h^' 



the PRISM thread description, in which the polymer 



chain is treated as an infinite thread of vanishing thickness. 



^0 



exp 



exp 



r 



(3) 



with = 3/{TTpP) with p the monomer site density. Moreover, is the length scale of 

density fluctuations defined as = + and = Rg/V^ is the length scale of 
the correlation hole.^ Eq-Q shows how the structure of the macromolecular liquid is fully 
specified once the monomer and the radius-of-gyration length scales are known. 

By assuming Gaussian intramolecular distributions for the monomer-monomer and 
monomer-com correlation functions, which is a well tested assumption, Yatsenko, et al. 
have shown that h"^{r) can be expressed in an analytical form in terms of molecular and 
thermodynamic parameters as^^ 



3^ 



Rn 
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Rn 



e"9 



(4) 



where erfc(a;) is the complementary error function and can be expressed in terms of the 
radius of gyration as Rg/{2'np*^j^ with p*^ = PchR^ being the reduced molecular number 
density. Although the PRISM-thread model employed at the monomer level description 
does not predict solvation shells, which are present in a polymeric liquid on the local scale, 
its description of the structure of the liquid at the center-of-mass length scale is correct. 
For example, Eq.Q recovers the isothermal compressibility of the melt at the k = limit^^ 
and it is thermodynamically consistent. This equation effectively maps the polymer melt 
into a liquid of soft colloidal particles of radius Rg and recovers, in its extension to the 
case of polymer mixtures, the formalism of Bhatia and Thornton for mixtures of colloidal 
particles.^^'^^'^*^ 
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When compared against simulation data, the analytical expression of Eq. Q reproduces 
well the com total correlation function obtained from UA-MD simulations for several poly- 
meric systems that have different degrees of polymerization, different monomeric structure, 
and different density and temperature.^^ Theoretical calculations of h'^'^ik) do not contain 
any adjustable parameter because the parameters entering Eq.Q are the ones used in the 
united atom simulation against which the formalism is tested. 

Starting from Eq. Q and enforcing for the closure equation the hypernetted chain ap- 
proximation, which is the most reliable closure for soft potentials, the effective intermolecular 
potential between the coms of a pair of soft colloidal particles, f'^^(r), is derived as^^"^'^ 

j3v'\r) = h'^r) - ln[l + h'\r)] - c'\r) , (5) 

with j3 = (ksT)'^ the inverse temperature of the system in Boltzman units. The direct pair 
correlation function, c'^^(r) is given in reciprocal space in terms of h'^'^{k) as 

Substitution of Eqjljand Eq|6]into Eq|5]gives the effective pair potential acting between two 
soft colloidal particles. This is the potential needed as an input to the mesoscopic simulation 
of the coarse-grained units and is the potential we use in our MS-MD. 

Because coarse-grained potentials result from the mapping of many-body interactions 
into pair interactions, through the averaging over microscopic degrees of freedom, they are 
parameter dependent. During the coarse-graining procedure, the potential acting between 
microscopic units, which is given by the Hamiltonian of the system, reduces to an effective 
potential, which is a free energy in the reference system of the microscopic coordinates. 
The coarse-grained potential so obtained contains contributions of entropic origin due to 
the microscopic, averaged-out degrees of freedom and is therefore state-dependent. This 
can be observed in the form of the total correlation function between coarse-grained sites, 
Eq.Q, which explicitly includes the structural parameters of the polymer, i.e. the radius-of- 
gyration and density screening length, as well as the thermodynamic parameter of the system 
number density. The temperature enters directly through Eq.(|5| and indirectly through the 
molecular parameters, such as Rg. 

The common procedure adopted to derive the effective potential between coarse-grained 
units is fully numerical. Once the level of coarse-graining is defined by the selection of the 
effective units, the pair distribution function between effective units, g{r), is numerically 
calculated. From g{r) the potential of mean force is derived as — ksT In g{r), and then 
used as an input to a coarse-grained simulation. By comparison against the UA simulation, 
corrections to the potential are obtained, and the new potential of mean-force is used as an 
input to a new simulation until full consistency is obtained between coarse-grained and UA 
simulations. Because numerically solved coarse-grained potentials have parameters that are 
defined through optimization against a specific UA simulation, they apply only for systems at 
thermodynamic conditions close to the one of the UA simulation, so they are not transferable. 
Moreover, because their form depends on the procedure used to derive them, if they are 
optimized against a given set of physical properties, it is not assured that they will reproduce 
other properties of the system. The numerical optimization of the potential against different 
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Figure 1: The total pair potential, which is input into our MS-MD simulation, is shown for 
different chain lengths of polyethylene (solid line). For comparison, the dashed line represents 
the corresponding mean-force potential. 



physical properties could lead to different effective potentials, numerically optimized, even 
for the same system. 

In this respect, an analytical potential such as the one we derive here has the advan- 
tage of being explicitly parameter dependent, and because it applies to systems defined by 
different values of the thermodynamic parameters, it is fully transferable. Moreover, because 
it is derived from a first principle approach, its properties are well defined. However, ana- 
lytical potentials, as the ones we derived here, can be obtained only for a limited number of 
systems, thus implying that in most cases the numerical optimization procedure is the only 
option in deriving the effective coarse-grained potential. In conclusion, both analytical and 
numerically evaluated potentials are useful as they answer complementary questions in the 
process of building optimized coarse-grained descriptions of macromolecular systems. 

Finally, it is important to notice that the mean-force potential is not the correct input 
potential for the mesoscale simulations, as the two agree only for systems at low density,^^ 
however it is a good first guess in the iterative numerical procedure commonly adopted. 
In Figure [T] we show for the systems investigated the direct potential, which is input to 
our simulation, and the mean-force potential. The two appear to be related, even if their 
numerical value and the extent of the attractive contribution is quite different. 
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3 Mesoscale Simulation 



Once the formalism to describe analytically the structure of a liquid of polymers in a coarse- 
grained representation is solved, and the effective potential between coarse-grained units is 
derived, it is possible to run directly mesoscale simulations of the coarse-grained systems. 
Our meso-scale simulations were performed by implementing classical MD simulations in 
the microcanonical [N, V, E) ensemble. We simulated liquids of point particles interacting 
through the soft potential resulting from Eqs.(|4])-([6]), which is repulsive at short range and 
slightly attractive at longer range. The interparticle potential differs depending on the 
chemical nature of the polymer investigated, as well as thermodynamic parameters of density 
and temperature. The systems we simulated included melts of polyethylene (PE) with 
= 44 (PE44), = 48 (PE48), = 66 (PE66), and = 96 (PE96) monomers. Other 
systems were melts of syndiotactic polypropylene (sPP), head-to-head polypropylene (hhPP), 
isotactic polypropylene (iPP), and polyisobutylene (PIB), all with 96 united atoms per chain. 
Each simulation was performed on a cubic box with periodic boundary conditions. We 
used reduced units in our simulation, such that all the units of length were scaled by Rg 
(r* = r /Rg) and energies were scaled by ksT. The density was fixed at the value reported in 
Table [Tj Temperature and radius-of-gyration were utilized for dimensionalizing the results 
obtained from the MS-MD simulations, after they were performed. 



Table 1: Simulation Parameters for Polyolefin Melts 



Polymer Number of CHx moieties T [K] p [sites/A'^] Rg [A] 



PE44 


44 


400 


0.0324 


10.50 


PE48 


48 


448 


0.0329 


10.88 


PE66 


66 


448 


0.0329 


13.32 


PE96 


96 


453 


0.0328 


16.78 


sPP 


96 


453 


0.0328 


13.93 


hhPP 


96 


453 


0.0336 


13.54 


iPP 


96 


453 


0.0328 


11.34 


PIB 


96 


453 


0.0357 


9.68 



The potential and its corresponding derivative entered as tabulated inputs to the pro- 
gram, and linear interpolation was used to determine function values not found in the sup- 
plied numerical data as the algorithm proceeds. In the initialization step, all particles were 
placed on a lattice, and each site was given an initial velocity pooled from a Mersenne Twister 
random number generator.''' Subsequently, the system was evolved using a velocity Verlet 
integrator. 

Equilibrium was induced in the ensemble by velocity rescaling. The desired tempera- 
ture in the ensemble was attained by a quencing procedure, where velocities were forceably 
rescaled at regular intervals during the equilibration stage. We checked that the system 
had reached equilibrium when we observed a Maxwell-Boltzmann distribution of velocities, 
a steady temperature, a stabilized Boltzmann if-theorem function, and a decayed transla- 
tional order parameter. At this stage, velocity rescaling was discontinued. We monitored 
the fluctuating steady temperature during the production stage to be sure that the sys- 
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tern remains in equilibrium. During the production stage trajectories were collected over a 
traversal of ~ 8Rg. 

A typical mesoscale simulation consisted of ~ 3000 particles, evolving for about 30, 000 
computational steps during ~ 4 hours, and performed on a single-CPU workstation. We 
stress that our mesoscopic simulations represent an underestimate in the computational time 
since these have not been subjected yet to a parallelization. The large-scale {Rg lengthscale) 
pair distribution function resulting from the mesocale simulation has been shown to be in 
quantitative agreement with the one obtained from a full UA simulation, so that no precision 
is lost on the large scale during the coarse-graining procedure. ^^"^'^ The mesoscale simulation, 
however, required a much smaller computational power than the full UA simulation. UA- 
MD simulations of polymer melts against which our MS-MD simulations are compared used 
~ 1, 600 chains with N — 96 units, where 1 ns of simulation required either 9 h of computer 
time on 512 processors of Sandia's ASCI Red Intel cluster or 25 h of computational time on 
64 processors on a Cplant DEC alpha cluster. ^^"^"^ The reduced computational time of the 
MS-MD allowed for improved precision in the large-scale regime as it was possible to increase 
the number of particles and the box size, with respect to the UA-MD, without dramatically 
affecting the computational time. 

In this paper we test our multiscale procedure by comparing the resulting monomer 
total distribution function against the UA-MD simulation just described. We call this simu- 
lations the "full UA-MD" . Those simulations were performed by Mondello et al.,^^ Jaramillo 
et al.,^^ and Heine et al.^''' 



Table 2: Comparison of Mesoscopic Simulations (MS-MD) with United Atom Simulations 
(UA-MD). Number of molecules ums in a box of length Lms in MS-MD; number of molecules 
riuA in a box of length Ljja in the full UA-MD 



System Ums Lms[^ 



riuA 



PE44 


4000 


175.7978 


100 


51.4033 


PE48 


4096 


181.6127 


400 


83.634 


PE66 


5324 


220.2708 


100 


58.5530 


PE96 


6859 


271.7286 


100 


167.270 


sPP 


4394 


234.2482 


1600 


167.27 


hhPP 


4096 


227.5623 


1600 


165.966 


iPP 


3375 


214.5239 


1600 


167.27 


PIB 


3375 


208.6244 


1600 


162.676 



To capture properties on the monomer scale, the multiscale procedure requires informa- 
tion on the local scale, which must be extracted from simulations with atomic level resolution, 
here UA-MD. However, the UA simulation input to a multiscale modeling procedure only 
requires a number of polymers of the order of , which is much smaller than the number 
of polymers simulated in the full UA-MD, as we will discuss further in the paper. Henceforth 
we refer to the UA-MD simulation input to the multiscale procedure as the "local UA-MD" . 
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4 Building the full structural description at the monomer 
scale using a multiscale approach 

In a multiscale procedure, simulations are performed on a hierarchy of models where the 
same system is coarse-grained at different levels of resolution. The success of the method 
depends on the extent to which the system properties are modulated around distinct length 
scales. In the description of a polymer melt, the two lengthscales that are involved are the 
monomer segment length and the polymer radius-of-gyration. As it is shown in the left 
panels of Figures [5]j8] for the polymers considered in this study the peak on the local scale 
and the global properties are well separated. The multiscale procedure would not work if the 
local and global features superimpose, which could happen for short chains of stiff polymers. 
For those systems, however, the full UA-MD simulations are not too time consuming and 
there is no need to apply coarse-graining or to adopt a multiscale procedure. 

In our mesoscale simulations polymers are represented as point particles interacting 
through a soft potential of the range of the polymer size. i.e. Rg. The MS-MD provide 
information for any lengthscale of the order and larger than the polymer radius-of-gyration. 
It should be noted that neither mesoscale nor UA simulations can accurately describe the 
system at a radius larger than half their simulation box length. However, because the box 
size can be considerably larger for mesoscale simulations than for UA-MD, the MS-MD can 
describe large-scale properties more efficiently than UA simulations (see Table [2]). 

On the local scale, however, mesoscale simulations do not provide any information. To 
recover the structure of the liquid at lengthscales smaller than or equal to Rg, it is necessary 
to combine MS-MD with UA-MD simulations. These UA-MD, however, only need to provide 
information on the local scale, and they are run on a small number of molecules, of the order 
of \/N, and are performed for a contained number of simulation steps. Even if UA-MD 
simulations have to be performed to provide the local information, the multiscale procedure 
efforts a net computational gain with respect to directly running a "full" UA-MD. 

On the global scale, small k, the monomer total distribution function is calculated 
from the data of the mesoscale simulation, using the inverse of Eq.g 



\k) 



h^^k) , (7) 



where mesoscale simulations provide h'^^{k) while UA MD simulations are used to determine 
the intramolecular form factors, uj'^^{k) and uj'^"^{k). 

Figure [2] and |3] show h"^"^{k) (solid line) calculated from the mesoscale simulation 
using Eq.([7| for several of the polymer melts of Table [l| where each chain has 96 CH^ units, 
with X = 1,2, or 3. The calculated h"^'^{k) (solid line) is compared to the full UA MD 
simulation (symbols). As expected, the theory compares well to UA MD simulations for 
the small k range (large r), and begins to diverge as k increases. This is a consequence of 
the coarse-graining procedure, Eq.([7]). Since monomers can superimpose with the center-of- 
mass, uj'^"^{k) approaches zero for large k, while u'^"^{k) has a finite value at contact due to 
the hard-core monomer-monomer repulsion. As a consequence the ratio: [uj"^"^ {k) / uj'^"^ {k)]'^ 
diverges at some large k value. 
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Figure 2: Plot of h"^"^{k) for different polymer melts. Solid line depicts h"^"^{k) calculated 
using Eq.([7]). Symbols represent data points from full UA MD simulation. 




Figure 3: Plot of h™'^{k) for different polymer melts. Solid line depicts h^^{k) calculated 
using Eq.([7]). Symbols represent data points from full UA MD simulation. 
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In general, h"^"^{k), calculated from Eq. ([T]), using data from mesoscale simulations and 
the intramolecular structure factors from UA simulations, accurately captures the structural 
properties of the polymer for small values of k (global structure), as it is shown in Figure 
[2j However, the same quality of agreement observed between mesoscale simulations and 
data of the monomer total correlation function for most of the cases is not observed in two 
specific polymer structures, namely for iPP and for PIB, see Figure |3} It is known that both 
PIB and iPP favor highly packed intramolecular structures because of the symmetry in the 
position of the sidechains in the monomeric unit.^^ Moreover, because of the presence of 
two side-groups in a PIB monomer, the number of monomers in a PIB chain containing 96 
united atoms, is equal to 24 monomers. With such a small number of monomers the PIB 
sample can hardly exhibit Gaussian statistics for the intramolecular distribution. Because 
the assumption of having a Gaussian statistics is the starting point in our approach, ^^"^^ this 
accounts for the discrepancies between the coarse grained values and UA MD data that we 
observe for these systems. 

5 Combining United Atoms and Coarse-Grained sim- 
ulations: the crossover lengthscale 

Since the mesoscale simulation only captures global properties, it is necessary to determine 
the length at which local, intramolecular effects become significant and cannot any longer be 
discarded. This defines the crossover lengthscale for the UA simulation at which data from 
UA-MD and MS-MD simulations are combined. 

The extent to which intramolecular effects remain important on large lengthscales 
depends on the length of the polymer and its flexibility. If the length of the polymer is 
constant, stiffer polymers span a larger volume and therefore have a higher (lower) number of 
inter- (intra) molecular contacts than their more flexible counterparts. Furthermore, longer 
chains of polymers with identical chemical structure have a higher statistical number of 
interpenetrating molecules and a higher (lower) number of inter- (intra) molecular contacts 
than shorter chains. The iPP and PIB samples presented here are characterized by efficient 
intramolecular packing, which is due to the particular monomeric structure: for these chains 
intramolecular interactions are dominant over intermolecular ones for an extended range of 
lengthscales. 

The extent of intramolecular packing is quantified by calculating the ratio between the 
number of intra and total site/site contacts at a fixed radial distance, r, 

f ir) = ^-^^^ (8) 
where Ns{r) is the number of intramolecular contact sites defined as 

Ns{r)=A7vp[ r'^uj'^'^{r')dr' . (9) 



The total number of sites in a given volume, Ntotai{f)i is given by 



Ntotai{r) = ^npr^ (10) 
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Figure 4: Fraction of intramolecular contacts, /^(r), for melts of polymers with different 
monomer architectures. 

Figure |4] shows fs{r) vs. r for different monomer architectures. As expected, iPP 
(dotted hne) and PIB (dot-dashed hne) show higher values of /s(r), whereas the stiffest 
molecule, PE (sohd line), exhibits the lowest value at a given radius. 

Following the reasoning presented at the end of this section, we assume that a fraction 
of intramolecular contacts equal to /s(r) = 0.025, which corresponds to 2.5% of the total 
sites involved in intramolecular contacts, defines the value of the crossover distance. This 
corresponds to a distance r that is unique for a given polymer chain, and it defines the 
corresponding value in reciprocal space, k = 27r/r, at which the two simulations have to be 
combined. 

Figures |5]j8] display the comparison of the total correlation function obtained with our 
procedure (Eq. ([T])) and data from extended UA simulations. The left panel of the Figures 
show the total correlation function in reciprocal space, h'^"^{k), for various polymer melts, 
(sPP, hhPP, iPP, and PIB) and for varying chain lengths of polyethylene (PE). /i™™(A;) is 
obtained by combining the large scale data (small k) from the mesoscale simulations with the 
small scale data (large k) from united atom simulations. The vertical dashed line in Figures 
[UIH] indicate the radius at which the MS-MD is combined with the UA MD simulation. Note 
that this occurs at an intermediate distance between the local peak and the global feature, 
so that neither information about the long range structure, nor the local structure is lost in 
extrapolating the connection between the two curves for h"^'^{k). 

The right panels of Figures [sjjs] show the total correlation function, /i™'"(r), for polymer 
melts, obtained by taking the Fourier transform of h"^"^{k). The function, /i™'™(r), provides 
a complete description of the liquid structure and thermodynamic properties. In performing 
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Figure 5: Left: Plot of the total correlation function, h"^"^{k), for polymer melts of (a) 
sPP and (c) hhPP, obtained by combining mesoscale and UA MD simulations. Mesoscale 
simulation depicts h^'^{k) over the small k range whereas UA simulation provides data over 
the large k range. The dashed line indicates the point at which the two simulations were 
combined. The inset depicts the local peak. Right: Plot of the related /i'""*(r), the total 
correlation function in real space for (b) sPP and (d) hhPP. The solid line depicts /i™'™(r) 
calculated using our multiscaling approach and the open circles represent data points from 
UA MD simulations. 
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k (A"') r (A) 

Figure 6: Same as in Figure [s] for polymer melts of (a) iPP and (c) PIB 

the discrete Fourier transform of the multiscale total correlation function, a sampling step of 
Ak = O.OIA"^ was used, which defines the resolution of the transformed function. Since we 
are combining two separate data sets, Ak must be large enough so that the Fourier transform 
is not affected by the small discontinuity at the point of intersection. As long as the interval 
of the discontinuity is of the same order as that of the sampling step there is no effect on 
the Fourier transform. 

The calculated values for /i'"™(r) using our multiscaling approach (solid line) are pre- 
sented along with data from the full UA MD simulations (symbols). The proposed method 
gives results in good agreement with the full UA simulation data and correctly captures all 
of the relevant structural features of the liquid, including solvation shells and the correlation 
hole observed in polymers. 

Once this procedure is adopted, the iPP and PIB samples also give a good agreement 
with the full UA-MD simulation. Because in these samples, and in PE44, the chains are 
short, the total correlation function still shows a small numerical error. The error is clearly 
displayed by the function in the r < d regime, with d the hard-core diameter of the repulsive 
inter-site interaction. This is the region of the excluded volume, and here it should be found 
that h{r) = -1. 

The local chain structure, at large-fc, is represented by the peak in the insets of Figures 
[5]j8| which has a shape that depends on the monomeric structure as well as on the thermo- 
dynamic parameters of density and temperature, as the local chain packing is affected by 
those quantities. Figure [9] shows superimposed the local peaks for polyethylene chains with 
increasing degree of polymerization at constant temperature of 448 K. Because the chains 
have the same monomeric structure, and the thermodynamic parameters for each sample are 
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Figure 7: Same as in Figure |5] for polymer melts of (a) PE44 and (c) PE48 



17 





g-20 
E 

% -40 
-60 





■I-'"' 

2 

■ ■ 

• 1 










■ 

-2 






■ 

i 


1 1.5 2 







1 

k (A ^ 





g-20 
I -40 
Q--60 
-80 

















;■ -2 

■ 


■ 




■ 

L 


1 1.5 2 : 


0.5 


1 1.5 


2 




k (A 



Figure 8: Same as in Figure IS] for polymer melts of (a) PE66 and (c) PE96 



close in value, the local peaks superimpose. This leads to a further computational gain for 
the multiscale procedure, with respect to performing the full UA-MD simulation. Because 
the local properties are largely independent on the global scale properties for samples with 
long polymer chains, the multiscale procedure allows one to obtain the local scale properties 
for all samples just from one local UA-MD, which can be performed on a melt of short 
polymer chains, e.g. N ^ 40. 

Once the pair distribution function is obtained from the multiscale procedure any ther- 
modynamic quantity of interest can be numerically derived, including the equation of state, 
the virial coefficient, internal energy, compressibility, and others. Clearly the degree of 
precision of the multiscale total distribution function will affect the precision of the result- 
ing thermodynamic quantities, so that care has to be taken in the choice of the crossover 
distance where local and global information is combined to optimize the agreement between 
UA-MD and multiscale data. 

The crossover distance, at which the UA-MD ceases to provide relevant information on 
the local scale, is determined from the chosen value of the fraction of intramolecular contacts. 
In all the calculations reported above we adopted the value /s('") = 0.025. To determine 
this value, we investigated the extent of agreement between the total correlation function 
calculated through the multiscale procedure and the one from the full UA-MD data as a 
function of fs- 

r) from UA simulation and h"^^(r) calculated 



Figure 10 shows a correlation plot of 



from our theory when local UA-MD data were combined with MS-MD data at a radius 
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Figure 9: Plot of the local peak for the total correlation function, h"^"^{k), for polymer melts 
of PE at increasing degree of polymerization, (A^ = 48, 66, 96) from local UA simulations at 
T=448 K. 

corresponding to fs{r) = 0.05 (top) and /s(r) = 0.025 (bottom). Shown are the correlation 
plots for the stiffest chain, PE, and for the two most flexible ones, PIB and iPP. The degree of 
correlation between theory and simulation increases as /s(r) decreases, as expected; however 
already at /^(r) = 0.025 there is a very high degree of positive correlation between the two. 
The correlation coefficients between theory and simulation for each polymer architecture are 
shown in Table [s] As /s(r) decreases, the mesoscale picture of the liquid becomes increasingly 
accurate. However, /^(r) only becomes appreciably small for a large radius, r, and so there 
is a tradeoff between computational time in performing UA MD simulations and obtaining 
an accurate picture of the liquid structure. 

Because the correlation coefficient is close to unity in all polymer systems tested for 
fs{r) = 0.025 and does not improve significantly for smaller values of fs{f), we propose that 
this value serves as a general rule of thumb for combining mesoscale and UA simulations. 

Table 3: Correlation coefficient for polyolefin melts 



polymer 


correlation coefficient 
at fs{r) = 0.025 


PE96 


0.9998 


sPP 


0.9999 


hhPP 


0.9998 


iPP 


0.9998 


PIB 


0.9995 
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Figure 10: Correlation coefficient for h"^"^{r) obtained from the full UA-MD simulation and 
its theoretical value calculated with /^(r) = 0.05. Bottom panels: same as in top panels, but 
with fsir) = 0.025. 

A possible different criterion of selection of the crossover lengthscale could be to defined 
the point in which the UA-MD curve intersect the MS-MD (see Figure |2|. This leads to 
values of the distance r that are slightly larger than the ones established above, and for this 
reason requires more demanding local UA-MD, without adding precision to the calculated 
pair distribution functions. Furthermore, as mentioned above, for PIB and iPP, the curves 
never superimpose making an exact criterion of this sort for the combination of the two 
simulations unclear. By using the method adopted above of calculating the fraction of 
intramolecular contacts to determine the relevant crossover length scale, an exact distance, 
at which to combine the information, is readily obtained, and this distance does not depend 
on explicitly knowing the full UA MD simulation. 



Table 4: Total Number of Sites (A^siies) and Number of Molecules (n) included in a spherical 
volume of radius r, evaluated for fs{r) = 0.025 and 0.05 



System N,ites{0.025) n(0.025) A^,ites(0.05) n(0.05) 



PE96 


2785 


29 


995 


10 


sPP 


3547 


37 


1412 


15 


hhPP 


3518 


37 


1511 


16 


iPP 


3931 


41 


1758 


18 


PIB 


3729 


39 


1749 


18 
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The value of = 0.025 is reached for a distance r that is polymer dependent. 

The number of polymer chains statistically comprised in the volume spanned by r defines 
an average number of chains, n, that we need to consider in a local UA-MD simulation to 
produce good statistical information. This number can be calculated from Eq. ( 10 ) using the 
crossover radius and dividing the number of particles so obtained by the number of CH^, with 
X = 0, 1, 2, 3, units contained in a single chain. The number of sites, A*"^^, and the number 
of polymer chains, n, contained in the volume defined by the given radius of intramolecular 
interactions, calculated for value of /s(r') = 0.025, as well as for /s(r) = 0.05, are presented in 
Table HI The number of chains n is also the statistical number of macromolecules that need 
to be considered in a local UA-MD simulation and is of the order of the number of chains 
interpenerating a "tagged" polymer, This number of chains is one order of magnitude 
smaller than the number of molecules commonly used in the full UA simulations of polymer 
melts. 



6 Discussion and Conclusions 

In this paper we presented a new multiscale approach to simulate the structure of polymer 
melts in a computationally efficient way. We first identified the relevant lengthscales in which 
information has to be collected. For a liquid of uncharged homopolymer chains there are two 
lengthscales, which characterize all the structural properties of the macromolecular liquid, 
namely the monomer and the molecular radius of gyration. 

Molecular Dynamics simulations are performed for each system in the two different 
coarse-grained representations. At the monomer level, the liquid is well represented through 
an United-Atom (UA) description, where each CHx unit, with x = 1,2, or 3, corresponds 
to an effective unit. At the molecular level, the coarse-graining is obtained through a soft- 
colloid representation, where each polymer molecule is depicted as a point particle centered 
at the polymer center-of-mass, and interacting with other particles through a soft repulsive 
potential of the range of Rg. 

The two coarse-grained descriptions are formally connected through an Ornstein-Zernike 
equation, solved under the condition that the coms are fictitious sites, while the united atoms 
are real sites. An analytical solution of the OZ equation in terms of the total correlation 
function between the two representations has been derived, and this equation is used to 
combine local and global scale information in the proposed multiscale modeling procedure. 

Information on large scale properties is collected from the MS-MD, soft-colloidal, simu- 
lation, which extends to large-scale and covers long-time regimes, as the dynamics is speeded 
up because of the elimination of the internal degrees of freedom. The large-scale information 
has been shown to be quantitatively correct. However, because local degrees of freedom are 
averaged out, MS-MD simulations do not provide any information on local-scale properties. 

To obtain a description of the macromolecular liquid at the monomer lengthscale, MS- 
MD simulations have to be combined with local UA-MD simulations. A crossover lengthscale 
between the two representation has to be defined. This is equivalent to establishing at which 
small lengthscale the coarse-grained description ceases to be valid. The local degrees of 
freedom are relevant up to a characteristic lengthscale, which is different depending on the 
chemical structure of the monomer, but which can be evaluated using short-time small-box 
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united atom simulations. Once this crossover lengthscale is determined, it is adopted as the 
point at which global lengthscale information from MS simulations has to be combined with 
the local scale description from UA simulations, in reciprocal space. The resulting monomer 
total distribution function is Fourier transformed from reciprocal to real space to give the 
complete description of the system, through the total correlation function, h{r). 

We show that the resulting function, tested against data obtained from a full UA- 
MD simulation, exhibits quantitative agreement. Furthermore, a comparison of the needed 
computational time and resources of the multiscale procedure, against the requirements of 
the full UA simulation for the same system, highlights several favorable points. 

In general, simulations of liquids necessitate the inclusion of a high number of molecules, 
n, to achieve a good statistical description, minimize finite-size effects, and approximate 
the thermodynamic limit (n — >■ oo and y — > oo while density, p — n/V, is constant). 
Moreover, macromolecules usually comprise a large number of monomers A^, so that the 
overall computational time to simulate polymeric liquids increases rapidly, scaling as {nN^, 
if we assume the conventional scaling of computational steps with the square of the number 
of particles. This feature poses some limitations on the total number of polymers, their 
length, or the longest timescale that can be simulated. 

Other techniques can be adopted to decrease the exponent in the scaling of the compu- 
tational time, but these procedures usually cannot be translated into parallel calculations. 
Even with the most powerful resources, long-time UA-MD simulations of polymer melts usu- 
ally can treat a maximum of 500 chains, for polymers comprised by less than 100 monomers. 
However, the best statistical description is reached for simulations of thousands of molecules, 
which is very time consuming and impossible to perform for long polymers. 

A multiscale procedure is advantageous, because the mesoscale simulation provides the 
information on the large spatial and long-time scale, so that microscopic UA simulations 
can be limited in the number of particles that are monitored and in the computational time 
during which data are recorded. 

In our procedure the number of molecules to be included in the UA simulation is of 
the order of the number of chains interpenetrating a single macromolecule (0(\/]V)), which 
is about one order of magnitude smaller than the number of particles needed if we perform 
a full UA simulation. For example, if we assume a scaling proportional to the square of the 
number of particles, the proposed procedure requires only = N soft colloidal particles to 
calculate h'^^{k) and only {nNY = number of units in the local UA simulation, which is a 
dramatic gain with respect to the UA scaling of {nN^ where n >> A^ for the full simulation. 

If we assume a different scahng with the number of particles, for example because of 
adopting a more advanced sampling rule, the multiscale approach is still more convenient 
in the mimber of computational steps than the full UA MD simulations, since the gain in 
computational time when comparing the two procedures is due to the change in the total 
number of molecules to be tracked in the UA simulation, i.e. from n ^ A" in the full UA 
simulation to n oc y/N in the local UA simulation. 

Finally, computational time is further reduced when a family of homologous chains that 
only differ because of the degree of polymerization are under study: in this case local UA-MD 
only have to be performed for the shortest chain, as discussed in the previous section. 

Having reduced the computational time necessary to obtain local information, the 
global scale description results from a mesoscale simulation, which can include a large number 
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of particles and large box size without dramatically increasing the computational time, thus 
improving the statistics of the system on the global scale. Table |2] shows the comparison 
between the total ensemble number (A^) and simulation box length (L) for UA simulations 
and MS simulations of different polyolefines. 

For example, the full UA simulation for polyethylene with 66 united atoms includes 
100 polymer chains, while our MS simulation includes over 5,000 polymers as soft colloids, 
providing a mesoscale liquid structure which is identical in the two simulations and a large 
scale description with improved statistics. This reinforces the point that only small UA 
simulations need to be performed to capture local properties, whereas the global properties 
may be deduced using a mescoscale approach which greatly extends the capabilities for 
computer simulations of macromolecules. The possibility of simulating systems on very large 
lengthscales while including a large number of particles can be useful in treating systems with 
diverging correlation lengths, such as mixtures approaching their spinodal decomposition. 

In conclusion. Figures [5]j8] show that our method of combining mesoscale simulations 
at the coarse grained level with small UA simulations accurately captures both the large and 
small scale structure of a polymeric melt while remaining computationally advantageous. 
The approach is analytical and uses no adjustable parameters, and is therefore applicable 
to a large number of different polymer systems of varying chain length, density, and with 
different local chemical structure. As an example in this paper we investigate polyethylene 
chains with increasing degree of polymerization, as well as polymer chains with different 
monomeric structure. The same multiscale procedure is used for all of them, providing 
efficient calculations of monomer static properties across the many lengthscales of interest. 
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